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ABSTRACT 

When a shock wave interacts with a contact discontinuity, there may appear a reflected 
rarefaction wave, a deflected contact discontinuity and a refracted supersonic shock. 
The numerical simulation of shock wave refraction at a plane contact discontinuity 
separating gases with different densities is performed. Euler equations describing inviscid 
compressible flow were discretized using the finite volume method on unstructured 
meshes and WENO schemes. Time integration was performed using a third-order Runge- 
Kutta method. The wave structure resulting from regular shock refraction is determined 
allowing its properties to be explored. In order to visualize and interpret the results of 
numerical calculations, a procedure for identifying and classifying gas-dynamic 
discontinuities was applied. The procedure employed dynamic consistency conditions and 
digital image processing methods to determine flow structure and its quantitative 
characteristics. The results of the numerical and experimental visualizations were 
compared (shadow patterns, schlieren images, interferograms). The results computed are 
in an agreement with the theoretical and experimental predictions of a regular refraction 
of a shock wave on an inclined contact discontinuity. 
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Introduction 

A shock wave process (SWP) is a process of reorganizing a gas-dynamic 
system with the parameters f into a system with the parameters f , where f and 

f are a set of gas-dynamic variables before and behind SWP. These sets include 

kinematic, thermodynamic and thermo-physical variables that define flow 
parameters. Specifically, examples of the variables include kinematic 
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parameters (velocity, acceleration), thermodynamic parameters (pressure, 
density, temperature, entropy, enthalpy), and thermo-physical parameters (heat 
capacity, adiabatic index, viscosity), which can change during SWP. 

The intensity of SWP is generally characterized by a ratio of static 
pressures J - pi/p behind and before SWP (pi is the static pressure behind SWP, 
and p is the static pressure before SWP). In some cases, this is denoted as 
J=p/p , where p is pressure behind SWP. The values J c > 1 describe flow 

compression, while the values J c < 1 describe flow expansion (rarefaction). 

Gas-dynamic discontinuity (GDD) represents a certain idealization of an 
area where parameters change discontinuously, replacing the area with surfaces 
where gas-dynamic variables unevenly change (Uskov, Bulat & Arkhipova, 
2014). GDDs in supersonic flows can be either zero-order or first order. 
Examples of zero-order GDDs include compression wave, shock wave, and 
sliding surface, wherein gas-dynamic flow parameters (pressure, total pressure, 
velocity) have discontinuities. First-order GDDs are also known as weak 
discontinuities (discontinuous characteristics, weak tangential discontinuities). 
In first-order GDDs, the first-order derivatives of gas-dynamic parameters have 
discontinuities. 

Waves are different from discontinuities because waves have finite width 
and fill up the area between the leading and trailing fronts where gas-dynamic 

variables change from the values /to/. All waves and discontinuities can be 

divided into two large groups (Bulat & Uskov, 2014a). The first group includes 
contact discontinuities (characterized by surface separating gases with different 
densities but equal pressures and flow velocities) and tangential discontinuities 
(characterized by sliding surface separating flows with equal pressures but 
different velocities). Gas does not flow through such surfaces. Therefore, these 
surfaces are not considered as SWPs. The second group includes normal waves 
and discontinuities (these should not be confused with straight). As gas does 
flow through normal waves and discontinuities, they are considered as SWPs. 

Waves and discontinuities can interact with each other. There are three 
types of wave interactions (. Uskov, Bulat & Arkhipova, 2014b) 

- GDDs (shock waves and contact discontinuities) crossing each other; 

- isentropic Riemann waves interacting with each other; 

- isentropic waves interacting with GDDs. 

In addition, waves and discontinuities can interact with solid surfaces. 

Contact and tangential discontinuities cannot cross each other. The process 
of shock wave or isentropic wave interaction with a contact or tangential 
discontinuity is called refraction. This process is accompanied by a wave 
refracting and reflecting at a contact (tangential) discontinuity. Other cases of 
wave interaction are called interference. Refraction of a moving shock wave at 
an inclined contact discontinuity is examined in this study. 

Shock wave refraction 

Several physical phenomena involving gas flows result in shock waves. 
Shock wave refraction problems attract increasing interest because of a wide 
range of possible practical applications, such as shock wave interaction with 
ocean surface and bubble media, underwater burst, and liquid combustion. Non- 
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uniform flow regions (such as high temperature areas, interfaces, dust clouds) in 
front of a shock wave produce qualitative changes in flow development. 
Examples of these phenomena include distortion of shock wave fronts, formation 
of new shock waves, high-pressure jets, and large-scale vortices, boundary layer 
separation; and jet accumulation (Jahn, 1956; Znamenskaya et al., 2011). 
Interaction of an oblique shock wave with an interface results in the 
development of complex unsteady shock wave configurations (Abd-El-Fattah & 
Henderson, 1978; Henderson, 1989; Henderson, 1991; Nouragliev et al., 2005). 
Shock waves moving through a medium containing gas bubbles or fluid globules 
result in distortion of shock wave fronts, shock wave accumulation, and multiple 
vortices. 

Extant research provided a theoretical analysis of the formed shock wave 
configurations using shock polars (Samtaney & Zabusky, 1993). Boundaries of 
areas with expansion wave refraction modes and areas with the refraction of a 
reflected shock wave under the changing obliquity angle of a contact 
discontinuity were determined. Previous research also discussed vorticity 
generation, vortex structure evolution, and mesh convergence (Samtaney & 
Pullin, 1996). Studies also provided an exact solution to the problem of plane 
shock wave refraction at a contact discontinuity (regular case) (Samtaney, 2003; 
Wheatley, Pullin & Samtaney, 2005). They also discussed the suppression of the 
Richtmyer-Meshkov instability when a magnetic field is imposed. Exact 
solutions to the problem were provided for both one-dimensional and two- 
dimensional cases of regular shock wave refraction at a contact discontinuity. 

Various shock wave refraction modes using experimental and numerical 
methods were examined (Zeng & Takayama, 1996). Numerical calculations were 
made on the basis of a TVD scheme of first-order accuracy and an approximate 
solution to the Riemann problem. The data obtained corresponded significantly 
well to the results of measurements and theoretical analysis. The only exception 
to this included high Mach numbers of an incident shock wave. 

Two-dimensional calculations based on the Godunov-type scheme of second- 
order accuracy were described (Yang, 1992). Discontinuity of density was 
aligned at 75 degrees to the horizontal, with the Mach shock wave number of 
1.2. An MUSCL scheme was used to calculate flow structure at different angles 
of discontinuity obliquity (Delmont, Keppens & van der Holst, 2009). Various 
finite difference schemes were also used in the calculations (Pember & 
Anderson, 2001). 

A numerical study of the Richtmyer-Meshkov instability during the 
interaction of a strong shock wave with a linear or sinusoidal discontinuity of 
density was conducted (Samtaney & Meiron, 1997). Research also provided a 
visualization of GDD using various approaches (Samtaney & Zabusky, 2000). 
The refraction of a shock wave when it interacts with a near wall layer of heated 
gas was examined (Banuti, Grabe & Hannemann, 2011). Studies also discussed 
a wide range of issues associated with the interaction of shock waves with 
contact discontinuities and the formation of the Richtmyer-Meshkov instability 
(Brouillette, 2002). 

The refraction of a spherical shock wave at an air-water interface was 
investigated (Dutta et al., 2004). The data computed were compared according to 
the shock-capturing scheme using the level-set method. This provides means for 
identifying a contact discontinuity and examining its evolution. Previous studies 
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discussed the application of discontinuity identification schemes with respect to 
the refraction of shock waves at interfaces (Liu, Khoo & Yeo, 2001) The method 
of characteristics was used to solve the problem at mesh nodes located near an 
interface. The developed approach was applied to examine shock wave 
interaction with a gas bubble. The level-set method was used to determine the 
location of a phase interface (Osher & Fedkiw, 2003). The level-set function 
complied with the transfer equation, and had a zero value at an interface or 
values opposite in sign at the areas of different phases. This made it possible to 
describe changes in the position of the interface. A solution method for problems 
with contact boundaries was also proposed (Burago & Kukudzhanov, 2005). 

This study conducts numerical simulation of plane shock wave refraction at 
a contact discontinuity. The discontinuity is initially linear and separates gases 
with different densities. A case of regular shock wave refraction, when a triple 
configuration of GDD was formed, was investigated. This consisted of an 
incident shock wave, a reflected shock wave, and a transmitted shock wave. 
Numerical calculation data processed in the form of shadow patterns, schlieren 
images, and interferograms were compared with optical flow visualization data. 

Flow structure 

Various refraction modes emerge at the moment of shock wave incidence on 
a contact discontinuity positioned at a certain angle to the horizontal and 
separating gases with different densities. These are characterized by a refracted 
wave front and the formation of either an expansion wave or a reflected wave. A 
number of these configurations were discussed in a previous study based on the 
data obtained from a physical experiment (Jahn, 1956). 

As a result of the interaction of the incident shock wave (I) with the contact 
discontinuity (C), the transmitted wave (T) and the reflected wave (R) are 
formed as shown in Fig. 1. These are either shock waves or expansion waves, 
depending on the intensity of the incident shock wave and density drop at the 
contact discontinuity. 



Figure 1 . Shock wave refraction at a contact discontinuity 
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The velocity of a transmitted shock wave depends on the ratio of densities 
at a contact discontinuity o = p2/pl. Hence, the velocity of a transmitted shock 
wave is either higher (slow—fast refraction, S/F) or lower (fast—slow refraction, 
F/S) than the velocity of an incident shock wave, depending on the ratio of 
surfaces (Abd-El-Fattah & Henderson, 1978; Henderson, Colella & Puckett, 
1991). Velocity qualifies as the acoustic impedance of media. Word order 
determines the direction of wave propagation. In the slow-fast (S/F, o > 1) case, 
the reflected wave is an expansion wave, while in the fast-slow (F/S, o < 1) case, 
it is a shock wave. Each group (S/F and F/S) can be divided into a number of 
additional configurations depending on the intensity of the incident shock wave 
(Nouragliev, 2005). Previous researches observed 8 different refraction modes in 
the S/F case and 4 different refraction modes in the F/S case by changing the 
experimental conditions (Abd-El-Fattah & Henderson, 1978). The conditions for 
the existence of each mode as well as the conditions for the transition from one 
mode to another were also examined. The stability of the configurations was 
examined (Fang, Wang & Yuan, 2011). 

With respect to certain parameters, a transition takes place from regular 
refraction (when three waves converge at a point) to non-regular refraction 
(when a Mach stem is formed). Extant research provides a detailed description 
and analysis of various modes of shock wave refraction at a contact discontinuity 
(Henderson, 1989; Nouragliev et al., 2005). 

Computational domain and boundary conditions 

The interaction of a plane shock wave moving from left to right and having 
a Mach number of M with a contact discontinuity oriented at the a angle to the 
direction of the flow is considered. The contact discontinuity separates gases 
with densities pi (to the left of the discontinuity) and p 2 (to the right of the 
discontinuity). 

Calculations were made in the area [0, 1] x [0, 5] on an unstructured mesh 
containing about 10 6 nodes, with the nodes accumulated near interface. Fig. 2 
shows a computational domain. In the calculations, the shock wave Mach 
number was assumed to be M = 2, whereas the contact discontinuity obliquity to 
the x-axis was a = 45°. The contact discontinuity separated two gases with the 
ratio of densities o = 3. Computational gas with the adiabatic exponent y = 1.667 
was used in the calculations. The calculations extended up to the time point 
t - 0.72 (in this time interval, a sound wave in an undisturbed gas ran across the 
width of a shock tube). 


C 




M 



y 



Figure 2. Computational domain 
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At the initial time moment, gases separated by a contact discontinuity were 
considered in a rest, and the shock wave front was located at xs = 0.1 (near the 
left boundary of the computational domain). The Rankine-Hugoniot relations 
were applied to the shock wave front. Non-reflecting boundary conditions were 
applied to the inlet and outlet boundaries of the computational domain. At the 
bottom boundary, symmetry conditions were applied (solid boundary). 

The numerical method 

The numerical model was constructed on the basis of the solution of the 
unsteady Euler equations for inviscid compressible gas (Volkov, 2008). 

The Euler equations were discretized by the finite volume method and an 
explicit WENO scheme of third-order accuracy on an unstructured mesh. 
Previous studies provide examples of the application of WENO schemes in the 
simulation of supersonic flows of inviscid compressible gas on unstructured 
meshes (Bulat & Volkov, 2015a; Bulat & Volkov, 2015b). Convective flows were 
calculated independently for each direction using HLLC (Harten—Lax—van Leer- 
Contact) Riemann solver. Time integration was performed using a third-order 
Runge-Kutta method. A perfect gas equation of state with a constant ratio of 
specific heat capacities was used. 

Discontinuity identification 

The process of gas-dynamic discontinuity identification (i.e., identification of 
their type and position) involves considerable time and does not lack in 
subjectivity. A method to accelerate the process and increase its accuracy was 
proposed (Vorozhtsov, 1990). This involved calculating the density gradient and 
its orientation in the center of each mesh cell. Points with a density gradient 
above the average across the area were considered as points of discontinuity. 
GDD were classified by examining the nearest points for discrete analogs of the 
conditions observed at strong discontinuities. Eventually, each point belonged to 
a certain type of discontinuity, namely shock wave (normal or oblique), 
tangential discontinuity, contact discontinuity, or compression wave 
(Vorozhtsov, 1990). 

Numerical differentiation causes noise interference in the resulting image. 
The order of accuracy of a finite difference scheme is decreased to the first-order 
near the discontinuity. This is usually required to preserve the monotonicity of 
the finite difference scheme. However, this also causes interference (Volkov, 
2008). A smoothing operation was performed to eliminate the noise interference 
(Schalkoff, 1988). 

Contact discontinuity identification involves the application of the level-set 
method. The level-set function complies with the transfer equation, and it has a 
zero value at an interface or values opposite in sign at the areas of different 
phases. This makes it possible to describe changes in the position of the 
interface. The condition at the contact discontinuity £( t) - 0 makes it possible to 
track its position in time (on the left and right side of the discontinuity, it is 
assumed that £ = ±1). 

Results and Discussions 
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In the calculations, the obliquity angle of a contact discontinuity to the 
horizontal is assumed to be a = 45°. With the given problem parameters, a 
regular refraction mode was implemented, and a reflected shock wave was 
formed R. Samtaney, 2003. A refraction mode wherein the reflected 
discontinuity was an expansion wave was examined P. Delmont, R. Keppens & 
B. van der Holst, 2009. The case when a = 90° corresponded to a one¬ 
dimensional solution to the problem of shock wave interaction with a contact 
discontinuity. 

Fig. 3 shows the density contours at the time moment t - 0.72. Gas density 
changes from 2.5 to 8.2. The density distribution shows the location of the 
reflected shock wave and contact discontinuity. 



Figure 3. Contours of density 


Fig. 4 shows the numerical schlieren at the time t - 0.72. The Sobel filter 
was used to discriminate discontinuities from the solution. Despite the fact that 
the Sobel method is less susceptible to numerical noise, in this particular 
problem, the use of the Roberts and Sobel operators resulted in virtually 
identical images (Schalkoff, 1988). 



Figure 4. Numerical schlieren 
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Fig. 5 shows the numerical shadow pattern, and Fig. 6 shows the numerical 
interferogram at the time moment t - 0.72. Increasing the number of lines 
makes it possible to distinguish the flow details. 



Figure 5. Numerical shadow pattern 



Figure 6. Numerical interferogram 


Fig. 7 shows flow field visualization using different variables (pressure 
Laplacian, velocity divergence, entropy Laplacian). The images shown in the 
fragments a and b are rather similar. Conversely, using the entropy Laplacian 
in flow visualization (fragment c) results in a loss (invisibility) of weak gas- 
dynamic discontinuities. 
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Figure 7. Flow visualization using pressure Laplacian (a), velocity divergence (b), entropy 
Laplacian (c) 


Fig. 8 shows shock wave structures. The lines 1, 3, 5, 6, and 7 correspond to 
shock waves, while the lines 2, 4, and 8 correspond to contact discontinuities. 
The points T1 and T2 are triple points, where two shock waves and a contact 
discontinuity cross. 














(5£) INTERNATIONAL JOURNAL OF ENVIRONMENTAL & SCIENCE EDUCATION 


9035 



Figure 8. Shock wave structure formed as a result of shock wave interaction with a contact 
discontinuity 


Fig. 9 shows contact discontinuity front positions determined using the 
density Laplacian (Ap = 0) and the level-set method (just a part of the 
computational domain is shown). Using various approaches, a significant 
difference in the form of determined contact discontinuity occurred only at the 
area of the maximum curvature. The shock wave front position depends on the 
approach applied for discontinuity identification. This dependence is 
significantly lower than that of a contact discontinuity. 



Figure 9. Contact discontinuity front positions determined using the density Laplacian (line 
1) and the level-set method (line 2) 
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Conclusion 

The interaction of a shock wave with an interface separating two gases 
produces a rich variety of flow phenomena even in the simplest configurations. 
Numerical simulation of plane shock wave interaction with an inclined contact 
discontinuity separating gases with different densities was performed. Results of 
the numerical visualization of shock wave flow patterns were compared with the 
data obtained from optical observations (including shadow patterns, schlieren 
images, and interferograms). Digital image processing methods were used to 
identify GDDs. The form of a contact discontinuity determined using the level- 
set method was compared with that obtained using a density Laplacian. 

The code developed could be used to investigate wave pattern in more 
complex cases including transition from slow-fast to fast-slow refraction. The 
experimental results show that after reflection from the top wall, the contact 
discontinuity becomes unstable due to Kelvin—Helmholtz instability, causing it 
to roll up and form a Richtmyer—Meshkov instability. These instabilities are 
subject of intensive research and will be considered in future works. The 
turbulent regime posed a significant challenge, and new experimental data are 
required to generalize the time evolution of the interface over a wide range of 
initial conditions. 
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